This analysis is based on 68 samples obtained from both datasets (diet-microbiota)
metadata <- read.csv("../data/Procrustes_metadata.csv", header = TRUE, check.names = FALSE,
row.names = 1)
diet <- read.csv("../data/Procrustes_Diet.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota <- read.csv("../data/Procrustes_Microbiota.csv", header = TRUE,
check.names = FALSE, row.names = 1)
#microbiota <- microbiota %>% filter(rowSums(across(where(is.numeric)))!=0) %>% t()# Distance matrix of bacterial microbiota
jaccard_mic_1 <- hillpair(data = microbiota, q = 1) # distance matrix at order q1## The estimated time for calculating the 2278 pairwise combinations is 4 seconds.
microbiota_jaccard_1 <- as.dist(jaccard_mic_1$S)
jaccard_mic_2 <- hillpair(data = microbiota, q = 2) # distance matrix at order q1## The estimated time for calculating the 2278 pairwise combinations is 3 seconds.
set.seed(123)
# Order q=1
mantel(diet_jaccard, microbiota_jaccard_1, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard, ydis = microbiota_jaccard_1, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.206
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0349 0.0473 0.0560 0.0687
## Permutation: free
## Number of permutations: 999
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard, ydis = microbiota_jaccard_2, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1948
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0387 0.0511 0.0610 0.0764
## Permutation: free
## Number of permutations: 999
# PCoA of microbiota q=1
microbiotaPCoA_1 <- as.data.frame(cmdscale(microbiota_jaccard_1))
plot(microbiotaPCoA_1)# PCoA of microbiota q=2
microbiotaPCoA_2 <- as.data.frame(cmdscale(microbiota_jaccard_2))
plot(microbiotaPCoA_2)procrust <- procrustes(X = dietPCoA, Y = microbiotaPCoA_1, scale=TRUE,
symmetric = TRUE)
pro_test <- protest(dietPCoA, microbiotaPCoA_1, permutations = 9999)
pro_test##
## Call:
## protest(X = dietPCoA, Y = microbiotaPCoA_1, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.8268
## Correlation in a symmetric Procrustes rotation: 0.4161
## Significance: 1e-04
##
## Permutation: free
## Number of permutations: 9999
eigen <- sqrt(procrust$svd$d)
percent_var <- signif(eigen/sum(eigen), 4)*100
beta_pro <- data.frame(procrust$X)
head(beta_pro)beta_pro$UserName <- rownames(beta_pro)
beta_pro$type <- "Diet"
beta_pro$seasons<-metadata[,5]
beta_pro$species<-metadata[,1]
head(beta_pro)trans_pro$UsarName <- rownames(trans_pro)
trans_pro$type <- "Microbiota"
trans_pro$seasons=metadata[,5]
trans_pro$species=metadata[,1]
head(trans_pro)colnames(trans_pro) <- colnames(beta_pro)
toplot <- data.frame(rbind(beta_pro,trans_pro))
pval <- signif(pro_test$signif, 1)
#col1 = rgb(250/255,60/255,60/255)
#col4 = rgb(0/255,200/255,200/255)
#col8 = rgb(160/255,0/255,200/255)
#col10 = rgb(0/255,160/255,255/255)
# "#999999","#E69F00","#56B4E9", "#009E73","#F0E442procrustes_new <- ggplot(toplot) +
geom_point(size = 1.8, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#999999","#E69F00","#56B4E9",
"#009E73")) +
theme_classic() +
scale_x_continuous(limits = c(-0.1, 0.3)) +
scale_y_continuous(limits = c(-0.2, 0.2)) +
geom_line(aes(x= V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
annotate("text", x = 0.2, y = -0.2,
label = paste0("p-value=", pval), size = 4) +
xlab(paste0("PCoA 1 (",percent_var[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var[2],"%)"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 3 rows containing missing values (`geom_point()`).
## Warning: Removed 3 rows containing missing values (`geom_line()`).
## Warning: Removed 3 rows containing missing values (`geom_point()`).
## Removed 3 rows containing missing values (`geom_line()`).
procrust <- procrustes(X = dietPCoA, Y = microbiotaPCoA_2, scale=TRUE,
symmetric = TRUE)
pro_test <- protest(dietPCoA, microbiotaPCoA_2, permutations = 9999)
pro_test##
## Call:
## protest(X = dietPCoA, Y = microbiotaPCoA_2, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.8177
## Correlation in a symmetric Procrustes rotation: 0.4269
## Significance: 1e-04
##
## Permutation: free
## Number of permutations: 9999
eigen <- sqrt(procrust$svd$d)
percent_var <- signif(eigen/sum(eigen), 4)*100
beta_pro <- data.frame(procrust$X)
head(beta_pro)beta_pro$UserName <- rownames(beta_pro)
beta_pro$type <- "Diet"
beta_pro$seasons<-metadata[,5]
beta_pro$species<-metadata[,1]
head(beta_pro)trans_pro$UsarName <- rownames(trans_pro)
trans_pro$type <- "Microbiota"
trans_pro$seasons=metadata[,5]
trans_pro$species=metadata[,1]
head(trans_pro)procrustes_new <- ggplot(toplot) +
geom_point(size = 1.8, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#999999","#E69F00","#56B4E9",
"#009E73")) +
theme_classic() +
scale_x_continuous(limits = c(-0.1, 0.3)) +
scale_y_continuous(limits = c(-0.2, 0.2)) +
geom_line(aes(x= V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
annotate("text", x = 0.2, y = -0.2,
label = paste0("p-value=", pval), size = 4) +
xlab(paste0("PCoA 1 (",percent_var[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var[2],"%)"))
print(procrustes_new)## Warning: Removed 3 rows containing missing values (`geom_point()`).
## Warning: Removed 3 rows containing missing values (`geom_line()`).
## Warning: Removed 3 rows containing missing values (`geom_point()`).
## Removed 3 rows containing missing values (`geom_line()`).
Procrsutes by lizard species were done at q=1 as similar patterns were observed in the general procrustes at q=2.
#loading data
aeneus_metadata <- subset(metadata, Species_Scel == "Sceloporus aeneus")
diet_aeneus <- read.csv("../data/Proc_Diet_aeneus.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_aeneus <- read.csv("../data/Proc_Mic_aeneus.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_a <- hillpair(data = microbiota_aeneus, q = 1)## The estimated time for calculating the 78 pairwise combinations is 2 seconds.
microbiota_jaccard_a <- as.dist(jaccard_mic_a$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_a <- vegdist(diet_aeneus, method = "jaccard")
diet_jaccard_a <- as.dist(diet_jaccard_a)
# Mantel test
set.seed(123)
mantel(diet_jaccard_a, microbiota_jaccard_a, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_a, ydis = microbiota_jaccard_a, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.3412
## Significance: 0.031
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.196 0.283 0.348 0.399
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_a <- procrustes(X = dietPCoA_a, Y = microbiotaPCoA_a, scale=TRUE,
symmetric = TRUE)
pro_test_a <- protest(dietPCoA_a, microbiotaPCoA_a, permutations = 9999)
pro_test_a##
## Call:
## protest(X = dietPCoA_a, Y = microbiotaPCoA_a, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.7374
## Correlation in a symmetric Procrustes rotation: 0.5125
## Significance: 0.0734
##
## Permutation: free
## Number of permutations: 9999
eigen_a <- sqrt(procrust_a$svd$d)
percent_var_a <- signif(eigen_a/sum(eigen_a), 4)*100
beta_pro_a <- data.frame(procrust_a$X)
head(beta_pro_a)beta_pro_a$UserName <- rownames(beta_pro_a)
beta_pro_a$type <- "Diet"
beta_pro_a$seasons <- aeneus_metadata[,5]
beta_pro_a$species <- aeneus_metadata[,1]
head(beta_pro_a)trans_pro_a$UsarName <- rownames(trans_pro_a)
trans_pro_a$type <- "Microbiota"
trans_pro_a$seasons <- aeneus_metadata[,5]
trans_pro_a$species <- aeneus_metadata[,1]
head(trans_pro_a)colnames(trans_pro_a) <- colnames(beta_pro_a)
toplot_a <- data.frame(rbind(beta_pro_a, trans_pro_a))
pval_a <- signif(pro_test_a$signif, 1)
# Graph diet and microbiota
procreustes_aeneus_new <- ggplot(toplot_a) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#999999")) +
theme_classic() +
scale_x_continuous(limits = c(-0.2, 0.4)) +
scale_y_continuous(limits = c(-0.2, 0.4)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.13, label = paste0("p-value=", pval_a),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_a[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_a[2],"%)"))
print(procreustes_aeneus_new)## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 row containing missing values (`geom_line()`).
#Loading data
bica_metadata <- subset(metadata, Species_Scel == "Sceloporus bicanthalis")
diet_bica <- read.csv("../data/Proc_Diet_bica.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_bica <- read.csv("../data/Proc_Mic_bica.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_b <- hillpair(data = microbiota_bica, q = 1)## The estimated time for calculating the 120 pairwise combinations is 2 seconds.
microbiota_jaccard_b <- as.dist(jaccard_mic_b$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_b <- vegdist(diet_bica, method = "jaccard")
diet_jaccard_b <- as.dist(diet_jaccard_b)
# Mantel test
set.seed(123)
mantel(diet_jaccard_b, microbiota_jaccard_b, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_b, ydis = microbiota_jaccard_b, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.3353
## Significance: 0.006
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.135 0.182 0.219 0.284
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_b <- procrustes(X = dietPCoA_b, Y = microbiotaPCoA_b, scale=TRUE,
symmetric = TRUE)
pro_test_b <- protest(dietPCoA_b, microbiotaPCoA_b, permutations = 9999)
pro_test_b##
## Call:
## protest(X = dietPCoA_b, Y = microbiotaPCoA_b, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.7286
## Correlation in a symmetric Procrustes rotation: 0.521
## Significance: 0.0226
##
## Permutation: free
## Number of permutations: 9999
eigen_b <- sqrt(procrust_b$svd$d)
percent_var_b <- signif(eigen_b/sum(eigen_b), 4)*100
beta_pro_b <- data.frame(procrust_b$X)
head(beta_pro_b)beta_pro_b$UserName <- rownames(beta_pro_b)
beta_pro_b$type <- "Diet"
beta_pro_b$seasons <- bica_metadata[,5]
beta_pro_b$species <- bica_metadata[,1]
head(beta_pro_b)trans_pro_b$UsarName <- rownames(trans_pro_b)
trans_pro_b$type <- "Microbiota"
trans_pro_b$seasons <- bica_metadata[,5]
trans_pro_b$species <- bica_metadata[,1]
head(trans_pro_b)colnames(trans_pro_b) <- colnames(beta_pro_b)
toplot_b <- data.frame(rbind(beta_pro_b, trans_pro_b))
pval_b <- signif(pro_test_b$signif, 1)
# graph diet and microbiota
procreustes_bica_new <- ggplot(toplot_b) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#E69F00")) +
theme_classic() +
#scale_x_continuous(limits = c(-0.25, 0.4)) +
#scale_y_continuous(limits = c(-0.2, 0.25)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_b),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_b[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_b[2],"%)"))
print(procreustes_bica_new)# Loading data
gram_metadata <- subset(metadata, Species_Scel == "Sceloporus grammicus")
diet_gram <- read.csv("../data/Proc_Diet_gram.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_gram <- read.csv("../data/Proc_Mic_gram.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_g <- hillpair(data = microbiota_gram, q = 1)## The estimated time for calculating the 253 pairwise combinations is 2 seconds.
microbiota_jaccard_g <- as.dist(jaccard_mic_g$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_g <- vegdist(diet_gram, method = "jaccard")
diet_jaccard_g <- as.dist(diet_jaccard_g)
# Mantel test
set.seed(123)
mantel(diet_jaccard_g, microbiota_jaccard_g, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_g, ydis = microbiota_jaccard_g, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.4457
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0905 0.1235 0.1568 0.1910
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_g <- procrustes(X = dietPCoA_g, Y = microbiotaPCoA_g, scale=TRUE,
symmetric = TRUE)
pro_test_g <- protest(dietPCoA_g, microbiotaPCoA_g, permutations = 9999)
pro_test_g##
## Call:
## protest(X = dietPCoA_g, Y = microbiotaPCoA_g, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.5908
## Correlation in a symmetric Procrustes rotation: 0.6397
## Significance: 1e-04
##
## Permutation: free
## Number of permutations: 9999
eigen_g <- sqrt(procrust_g$svd$d)
percent_var_g <- signif(eigen_g/sum(eigen_g), 4)*100
beta_pro_g <- data.frame(procrust_g$X)
head(beta_pro_g)beta_pro_g$UserName <- rownames(beta_pro_g)
beta_pro_g$type <- "Diet"
beta_pro_g$seasons <- gram_metadata[,5]
beta_pro_g$species <- gram_metadata[,1]
head(beta_pro_g)trans_pro_g$UsarName <- rownames(trans_pro_g)
trans_pro_g$type <- "Microbiota"
trans_pro_g$seasons <- gram_metadata[,5]
trans_pro_g$species <- gram_metadata[,1]
head(trans_pro_g)colnames(trans_pro_g) <- colnames(beta_pro_g)
toplot_g <- data.frame(rbind(beta_pro_g, trans_pro_g))
pval_g <- signif(pro_test_g$signif, 1)
# Graph diet and microbiota
procreustes_gram_new <- ggplot(toplot_g) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#56B4E9")) +
theme_classic() +
#scale_x_continuous(limits = c(-0.25, 0.4)) +
#scale_y_continuous(limits = c(-0.2, 0.25)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_b),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_g[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_g[2],"%)"))
print(procreustes_gram_new)spi_metadata <- subset(metadata, Species_Scel == "Sceloporus spinosus")
diet_spi <- read.csv("../data/Proc_Diet_spi.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_spi <- read.csv("../data/Proc_Mic_spi.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_s <- hillpair(data = microbiota_spi, q = 1)## The estimated time for calculating the 120 pairwise combinations is 2 seconds.
microbiota_jaccard_s <- as.dist(jaccard_mic_s$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_s <- vegdist(diet_spi, method = "jaccard")
diet_jaccard_s <- as.dist(diet_jaccard_s)
# Mantel test
set.seed(123)
mantel(diet_jaccard_s, microbiota_jaccard_s, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_s, ydis = microbiota_jaccard_s, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1073
## Significance: 0.151
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.141 0.190 0.224 0.267
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_s <- procrustes(X = dietPCoA_s, Y = microbiotaPCoA_s, scale=TRUE,
symmetric = TRUE)
pro_test_s <- protest(dietPCoA_s, microbiotaPCoA_s, permutations = 9999)
pro_test_s##
## Call:
## protest(X = dietPCoA_s, Y = microbiotaPCoA_s, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.797
## Correlation in a symmetric Procrustes rotation: 0.4505
## Significance: 0.0789
##
## Permutation: free
## Number of permutations: 9999
eigen_s <- sqrt(procrust_s$svd$d)
percent_var_s <- signif(eigen_s/sum(eigen_s), 4)*100
beta_pro_s <- data.frame(procrust_s$X)
head(beta_pro_s)beta_pro_s$UserName <- rownames(beta_pro_s)
beta_pro_s$type <- "Diet"
beta_pro_s$seasons <- spi_metadata[,5]
beta_pro_s$species <- spi_metadata[,1]
head(beta_pro_s)trans_pro_s$UsarName <- rownames(trans_pro_s)
trans_pro_s$type <- "Microbiota"
trans_pro_s$seasons <- spi_metadata[,5]
trans_pro_s$species <- spi_metadata[,1]
head(trans_pro_s)colnames(trans_pro_s) <- colnames(beta_pro_s)
toplot_s <- data.frame(rbind(beta_pro_s, trans_pro_s))
pval_s <- signif(pro_test_s$signif, 1)
# Graph diet and microbiota
procreustes_spi_new <- ggplot(toplot_s) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#40A175")) +
theme_classic() +
#scale_x_continuous(limits = c(-0.25, 0.4)) +
#scale_y_continuous(limits = c(-0.2, 0.25)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_s),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_s[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_s[2],"%)"))
print(procreustes_spi_new)tree.diet <- read.tree("../data/tree.nwk")
comm<-(diet)
IDs= colnames(comm) # para obtener los nombres de las sp
tree.diet= ape::rtree(n=ncol(diet), tip.label = paste0(IDs)) # tener el árbol enraizado, con los nombres de los sitios y los nombres de las sp.
hill_phylo_parti_pairwise(comm,tree.diet, q = 0, show_warning = F) # Loading phylogenetic tree
phy_tree <- read_tree("../data/tree.nwk")
# Creating phyloseq object
SAM <- sample_data(metadata)
OTU <- otu_table(diet, taxa_are_rows=F)
phylo_physeq = phyloseq(OTU, SAM, phy_tree)
#Calculating phylogenetic distance
diet_dist_phy = distance(phylo_physeq, method = "uunifrac") %>% as.dist